High-order correlation effects in the two-dimensional Hubbard model 
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The electronic states of the two-dimensional Hubbard model are investigated by means of a 4-pole 
approximation within the Composite Operator Method. In addition to the conventional Hubbard 
operators, we consider other two operators which come from the hierarchy of the equations of motion 
and carry information regarding nearest-neighbor spin and charge configurations. By means of this 
operatorial basis, we can study the physics related to the energy scale of J = 4t 2 /U in addition 
to the one of U. Present results show relevant physical features, well beyond those previously 
obtained by means of a 2-pole approximation, such as a four-band structure with shadow bands and 
a quasi-particle peak at the Fermi level. The Fermi level stays pinned to the band flatness located 
at (7r,0)-point within a wide range of hole-doping (0 < 5 < 0.15). A comprehensive analysis of 
double occupancy, internal energy, specific heat and entropy features have been also performed. All 
reported results are in excellent agreement with the data of numerical simulations. 
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I. INTRODUCTION 

The discovery of high-T c superconductors promoted a 
largely diffused revival of interest in strongly correlated 
electron systems and fostered the study of many other 
transition-metal oxidesi Since the very beginning, the 
Hubbard model^ has received particular attention as it is 
retained a prototype for strongly correlated electron sys- 
tems and a minimal model to describe transition-metal 
oxides. In spite of the deceiving simplicity of its Hamil- 
tonian, a deep comprehension of all its physical features 
is still missing. In particular, it is very difficult to prop- 
erly describe the competition between kinetic, diagonal 
in momentum space, and potential, diagonal in direct 
space, energy terms. The main gross feature of the model 
is the splitting of the electronic band into two subbands 
divided by a gap of the order U. Nowadays, it is well- 
established that this feature can be understood in terms 
of well-known Hubbard operators within a 2-pole approx- 
imation. However, within this framework, inter-site cor- 
relations are poorly taken into account while they are 
universally recognized as essential to describe low-energy 
physics near half filling. For instance, Kampf and Schri- 
effer pointed out that, in highly correlated electron sys- 
tems, antiferromagnetic spin fluctuations play a funda- 
mental role in understanding features as pseudogap and 
shadow bands Recently, the remarkable progress in ex- 
perimental techniques has made possible to reveal such 
low-energy features in high-T c superconductors^ Anti- 
ferromagnetic fluctuations are often included within self- 
energy in a phenomenological way. Therefore, it is quite 
natural to wonder if, being these features inherent to 
the model Hamiltonian, it is possible to derive them 
microscopically by means of a suitable analytical tool. 
Numerical-simulation techniques, as Exact Diagonaliza- 
tion and quantum Monte Carlo (QMC), positively answer 
to this question^SiLSiS although they are applicable only 
to small clusters due to the exponential increase of the 
Hilbert space with the system size. Within the frame- 
work of projection methods, in order to describe inter- 



site correlations and related spin and charge fluctuations, 
it is necessary to take into account high-order operators 
which carry information regarding nearest-neighbor cor- 
relations. We will move along this way. 

In this paper, we investigate the electronic states of 
the Hubbard model by means of the Composite Op- 
erator Method (COM) 10 ! 11 ! 12 ! 13 ! 14 which has shown to 
be capable of describing the physics of many strongly 
correlated systems. On one hand, we can find many 
similarities between COM and the projection opera- 
tor method 15 i 16 i 17 i 18 i 19 i 20 i 21 i 22 and the spectral density 
approach (self-consistent moment approach) ,2i2£i2£i 2 £i2£ 
On the other hand, COM differs from these meth- 
ods as regards the conscious exploitation of the pres- 
ence of unknown parameters in the theory in order 
to put constraints on the representation where the 
Green's functions are realized. The Hubbard model has 
been widely analyzed by means of COM within the 2- 
pole approximation^i 2 *^ In this approximation, COM 
reaches a global agreement with numerical simulations 
regarding local and thermodynamic quantities. In order 
to go beyond the 2-pole approximation, it is necessary 
either to evaluate the dynamical corrections or to intro- 
duce high-order operators in the basis. As regards the 
former case, a fully momentum and frequency dependent 
self-energy has been evaluated by means of: two-site level 
operators within the non-crossing approximation, 28 29-30 
loop decoupling within both the self-consistent Born 



approximatio! 



.31,32,33 



and the iterative perturbation the- 



ory within the dynamical mean field theory^ As regards 
the latter case, Dorneich et al. have introduced in the op- 
eratorial basis, in addition to the conventional Hubbard 
operators, two nonlocal composite fields which describe 
the local electronic transitions dressed by the nearest- 
neighbor spin and charge fluctuations,— They have man- 
aged to reproduce the well-known four band structure in 
good agreement with the QMC data, but only at half fill- 
ing to which their formulation is unfortunately restricted. 
According to all this, we have decided to examine the 
model by means of a 4-pole approximation within the 
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Composite Operator Method. In addition to the conven- 
tional Hubbard operators, we have also considered in the 
basis other two operators coming from the hierarchy of 
the equations of motion. Within our formulation, we can 
evaluate the evolution of each subband up to the second 
order in the equations of motion and we have managed 
to perform the analysis with finite doping, which was 
out of the scope in Ref. 35. Our results present: a four- 
band structure, a quasi-particle peak at the Fermi level, 
shadow bands, band flatness at (ir, 0)-point, a Fermi level 
located around the (ir, 0)-point and robust with respect 
to the hole doping (5 < 0.15). A comprehensive study, 
and comparisons with numerical simulations present in 
the literature, has been performed as regards: density of 
states, band dispersion, double occupancy, internal en- 
ergy, specific heat and entropy. 

The manuscript is organized as follows. In the next 
section (Sec.[Hj|, we nx the notation regarding the Hamil- 
tonian and give the general framework of the Compos- 
ite Operator Method. In Sec. IIIII the density of states 
and the dispersion relation are computed and discussed, 
also in comparison with some quantum Monte Carlo re- 
sults. A detailed comparison with numerical simulations 
for many local and thermodynamic quantities is given in 
Sec. IIVI Section contains summary and conclusions. 
Some detailed derivations of formulas used in Sec. IIIII are 
given in the appendix. 



II. MODEL AND FORMULATION 

The two-dimensional grand-canonical Hubbard Hamil- 
tonian H reads as 

H = H-fiJ2 n ^ (!) 
H = Y,t ii 4(iWU) + Uj2Mi)H(i) (2) 

ijCT i 

where c\(i) and c a (i) are creation and annihilation opera- 
tors, respectively, of electrons with spin a at the site i [i = 
(i, t)]. n a (i) — cJ.(i)co-(i). [i is the chemical potential. 
<ij = — 4tojj. a(k) = .F[ay] = | {cos(fc x a) + cos(fcj,a)}. 
a is the lattice constant. T is the Fourier transform. 
U is the on-site Coulomb repulsion. Here, we consider 
nearest-neighbor hopping only. We define the following 
operatorial basis 



1pBcr{i) 
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(i)J 
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Vso- 



(3) 



where £ ff (£) = (1 - n_ CT (i)) c a (i) and r) a (i) = n- a (i)c a (i) 
are the usual Hubbard operators that describe the tran- 
sitions of the local electron number n(i) — <-> 1 and 
1 <-> 2, respectively. The equations of motion of the com- 



ponents of ipAii) read as 
. d 



i^&(0 = -A*&(i) - 4*K(i) + 7r ff (i)} (4) 
d 



at 



(•5) 



where 



n a (i) = -n- a (i)c%{i) + cl a (i)c a (i)c°L a {i) 

+ c (7 (i)c" t .(i)c_ (r (£) (6) 

C^)=£>yC CT (j,i). (7) 
j 

Now, we can divide ir a (i) into two operators, ~K„(i) — 
Cscr(*) + ??so-(i), similarly to what we have done with c a {i), 
c a (i) — ^(i) + r] fT (i). That is, we choose the components 
of ipsii), Cscr(*) and ri sa (i), among the eigenoperators of 
the two-site Hubbard model^ and of the local interaction 
term of the Hamiltonian QJ. Then, they read as 

= -n- a {i)g{i) 

+ cl a (i)c a (i)CAi) + c^)7f_i(i)c^(i) (8) 

+ cU*M*)^M + Ccr (i)CUi)c- a (i). (9) 

It is clear now that ipB(i) describes nearest- neighbor com- 
posite excitations and carries information regarding sur- 
rounding spin and charge configurations^ According to 
the way we have chosen them^& £<r(£) and £,sa(i) belong 
to the energy class of the lower Hubbard subband and 
rj a (i) and rj sa (i) belong to the energy class of the upper 
Hubbard subband (see Eqs. (£), ©, JSUJ) and |5T}). 

Within the Composite Operator Method, once we 
choose a n-component operatorial basis ip, its equation 
of motion reads as 



d 



Y,e(i,j)il>Q,t)+6j(i) 



(10) 



where e(k) = m(k)7 1 (k) is a n x n matrix with 

I(k)=T({iP(i,t),^(j,t)}) (11) 

(12) 



d 



m(k) = F({if^(i,t),^(j,t)}}. 

The spinor notation is understood and (• • • ) denotes 
the thermal average in the grand-canonical ensem- 
ble. The expression of e(k) comes from the request 
({Sj(i, t), ip'(j, t) }) = 0. This constraint assures that the 
residual term Sj(i) contains only the dynamical correc- 
tions in terms of orthogonal components to the chosen 
basis. Neglecting Sj(i) gives the n-pole approximation 
for the retarded electronic Green's function G(u>,k) = 
T(Rm)^(j)]): 



G(w,k) 



E 



ffi(k) 



(13) 
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where Ei(k) are the eigenvalues of e(k) and cr,(k), the 
spectral functions, can be computed in terms of the eigen- 
vectors of e(k) and of the elements of I(k)M^ 

In the last years, the 2-pole approximation, that is, 
tp = i/iA) has been analyzed in great detail^ In the 
present paper, we perform a 4-pole analysis by enlarg- 
ing the operatorial basis with the introduction of tj) B - In 
this case, 7(k) reads as 



/(k) = 



/ Iaa (k) I AB (k) 

\l AB (k) Ibb (k) 

( hi 

7 22 



/i 3 (k) 
7 24 (k) 



7 X3 (k) \ 
7 24 (k) 



/ 33 (k) 7 34 (k) 
7 34 (k) 744 (k) / 



where 



with 



7llcr = 1 - (n- a ) 

I 13a (k) =A a - ar(k) ((n_ CT ) - Pa ) 
7 24(T (k) = -A CT - a^k)^ 



(14) 

(15) 
(16) 
(17) 
(18) 

(19) 
(20) 



The detailed expressions of the elements in the block 
IbbO*.) are rather complicated and reported in Ap- 
pendix. In order to effectively perform calculations, we 
have decoupled these elements by paying attention to 
the particle-hole symmetry enjoined by the Hamiltonian. 
Under this transformation, we have 



- (-iMM 

^(i)^(-i)^t(i ) + (-i)^t CT (i) 
-(-1)^(0 + (-i)U(O- 



(21) 
(22) 
(23) 
(24) 



After these relations, we have the following constraints 
on the 7o-(i,j) = ({^(i, t), V4-(ji *)}) matrix elements 



7n(i,j) - 


K-i) 


+j 7 22 (j,i) 




(25) 


722 (i,j) - 


H-i) 


+j 7n(j,i) 




(26) 


733(i, j) - 


*(-i) 


+j {/ 22 (j Q ,i a ) 


+ 7 24 (r,i)} 








+J{7 42 (j,i a ) + 


7 44 (j,i)} 


(27) 


^34(1, j) - 


*(-i) 


+j 7 34 (j,i) 




(28) 


7i4(i, j) - 


,(-1) 


+j {/n(r,n 


+7 13 (r,i)} 




- 


-(-1) 


+j U3i(j,i Q ) + 


7 33 (j,i)} 


(29) 



where i a stands for the nearest-neighbor sites of i (e.g., 
7i 3 (j",i) = ctji/i 3 (l, i)). It is worth noticing that our 
decoupling procedure exactly satisfies these relations. 



Now, we introduce a new operator 

■tpBa(i) = 

with 




UW-A 3 i(-iv)^W 

r) sa (i) - A i2 {-'^)ria{i) 



^[A 3 i(-iV)] = A 31 (k) = ^|^ 

^[^ 42 (-iV)] = J 4 42 (k) = ^iM. 

J 22 

This operator gives 7 in a block-diagonal form 

\ 



(30) 

(31) 
(32) 



= ( lAA(k) 

I Ibb (k) 



( hi 

7 22 





V 



7 33 (k) 7 34 (k) 
7 34 (k) 7 44 (k) J 



(33) 



where 



7 33 (k)= 7 33 (k) 
7 34 (k) = 7 34 (k) 
7 44 (k) = 7 44 (k) 



A 2 3 (k) 



'11 



^ 2 4 (k) 



1-22 



(34) 
(35) 
(36) 



Hereafter, we will use this new operator ^Baii) instead of 
^Ba{i) as it allows to more easily distinguish the contri- 
butions to single-particle properties of type B operators 
from those of type A. Then 



Vv(«) 



ipAv(i) 



(37) 



By means of this new basis, we have the matrix e 
to7 _1 , 



e(k) 



caaQs.) CAs(k) 
es^k) £Bs(k) 

( en (k) e 12 (k) 
£21 (k) e 22 (k) 



£31 (k) e 32 (k) 
\ 641 (k) e 42 (k) 



eia(k) ei 4 (k)\ 
£ 23 (k) e 24 (k) 



£33 (k) £3 4 (k) 
£ 43 (k) £ 44 (k) / 



(38) 



where 



e u (k) = -/*-4t(a(k) 



As 00 



'11 



(39) 




(co-juj/t (0)~/l)/t 




(0, 0) 



(71, 0) k (JC, 7C) 



(0, 0) 



FIG. 1: The density of states and the corresponding dis- 
persion relation (solid line) at U/t = 8, n = 0.90 and 
T/t = 0.01. The width of the dispersion line represents the 
intensity of peak. The 2-pole solution (dashed line) and the 
non-interacting case (dotted line) are also reported. 



ei2 (k) 
621 (k) = At 



At I a(k) 
A 3 (k) 



I 22 



e22 (k) = - Ai + C/ + 4<^iM 

-*22 

ei 3 (k) = eu(k) = -£23(k) = -£24 (k) 

e3l(k) = _ u !sn±MiL 
eM(k) = 4t l2ffl±liffl 

^22 



= -At 



e 4 l (k) = -4t : 



^ 43 (k)+I 44 (k) 

hi 

£42(k) = 4i ^w+M 
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(42) 
(43) 
(44) 
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(46) 

(47) 



and 
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FIG. 2: Doping dependence of density of states and dispersion 
relation at n — 0.88, 0.90 and 0.92. Other external parameters 
are the same as in Fig. Q 



where 



d - 

m BB (k) = F {{i-^B<j{i-,t),il} J{ B(J (],t)}) 

_ f ra 33 (k) m 34 (k) \ 
\ m 34 (k) m 44 (k) I 



(49) 



e BB (k) = m BB (k)/ B ]j(k) 



(48) 



The operator ipBtr(i) provides two-site excitations and 
represents the natural extension of %pA a {i), which de- 
scribes one-site excitations, in the sense of a series expan- 
sion over finite-cluster excitations. Equations of motion 
of ipBa(i) are much more lengthy and have a much more 
complex form with respect to those of ipAa-(i) as they 
contain many three-site composite operators. The appli- 
cation of a systematic projection/truncation procedure, 
as the one applied to the equations of motion of tpAa-(i)t is 
just unfeasible in this case as, besides to be very lengthy, 
it would lead to the appearance of a plenty of unknown 
correlation functions in the energy matrix. In order to fix 
all these correlation functions, we would be forced to use 
some decoupling and would completely lose any possibil- 
ity to control the approximation. Then, we have opted 
for a controlled, at least in philosophy, approximation 
at the level of equations of motion and decided to ne- 
glect irreducible three-site operators by paying attention 
to evaluate exactly all one- and two-site components^!. 
This choice has only one obvious drawback: the two-site 
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FIG. 3: The dispersion relation at U/t = 8, T/t = 0.5 and FIG. 4: Double occupancy (solid line) versus U/t at (a) n = 
n = 0.75 and 0.94. The 2-pole solution (dashed line) and 8/9 and T/t = (Lanczos data (circle) from Ref. H| and (b) 
QMC data (circle) of Ref. Q are also reported. n - 1.0 and T/t = 1/16 (QMC data (circle) from Ref. 



correlations, not damped by three-site processes, result 
quite enhanced. 

(50) 

i|^.<rW - (-i" + ^)^« + 4i|^(i)+CW} • 

(51) 

Equations of motion 1|50[) and (|51|l give a simplified 
form of ess(k) 

e 33 (k) = -,i + 4i (a(k) + (52) 

e 34 (k) = At ^2a(k) + (53) 

e 43 (k)=4t( a( k)-^)) (54) 

e 44 (k) =-/x + c7-4t%^. (55) 

^22 

The correlation functions appearing in / and e, except 
for p a , can be now self-consistently determined through 



the retarded Green's function 

W f > = (J^f J J dkdu [1-Mu)] (-^[G(w,-k)]j 

(56) 

where Jf{u) is the Fermi distribution function. 

The parameters p a are out of the scheme of the present 
formulation as they cannot be directly connected to the 
Green's function under study. They describe nearest- 
neighbor spin, charge and pair correlations and, accord- 
ing to their actual values, play a fundamental role in 
determining the behavior of the system (antiferromag- 
netic character of the band dispersion, presence of metal- 
insulator transition, etc.)^ The use of operators not sat- 
isfying canonical commutation relations leads to the ap- 
pearance of unknown correlation functions in the formu- 
lation. The presence of these unknown correlation func- 
tion should not be seen as an inconvenience, as many 
other formulations do, but as an opportunity given by 
the method to implement exact relations dictated by 
symmetries and/or general principles which are not au- 
tomatically satisfied, that is, which are no more em- 
bedded in the Hilbert space of the composite operators 
whose Green's functions are computediiiii According 
to this, we will evaluate p a by means of the algebra 
constraint»i2ii4 {Qr} a ) = 0. These constraints ensure 
that no state referring to a irip/e-occupied site (obviously 
forbidden by the Pauli principle) is taken into account in 
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FIG. 5: Internal energy and specific heat at U/t = 4 and 
n = 0.75 and 0.90 (solid line). Finite temperature Lanczos 
data (dashed line) from Ref. |4(J and QMC data (circle) from 
Ref. |4l| are also reported. 



the averaging procedure and give the possibility to make 
the correct spin/particle counting. This allows to fulfil 
the particle-hole symmetry and to correctly describe the 
virtual processes at the basis of the low-energy processes 
(J scale of energy). A comprehensive analysis by means 
of the two-pole approximation has shown that the use of 
algebraic constraints provides very good agreement with 
the numerical simulation results well beyond the conven- 
tional Hubbard-I and Roth's decoupling schemed 



III. BAND STRUCTURE 

Figure ^ shows the density of states and the corre- 
sponding dispersion relation for U/t = 8, n — 0.9 and 
T/t = 0.01. As a first consequence of choosing a ba- 
sis constituted of four operators, we have a four-band 
structure. Together with the usual Hubbard splitting of 
the non-interacting band with the appearance of a gap 
of the order U between the Hubbard subbands, we can 
clearly observe other two dispersion lines in the lower and 
higher energy regions (w/t ~ —4 and 8). They can be 
interpreted as shadow bands coming from the antiferro- 
magnetic nature of the spin fluctuations. In fact, they 
show a tendency towards a doubling of the zone through 
a mirroring of the original dispersion lines. Another re- 
markable feature is the presence of a well developed peak 
structure at the Fermi level which comes from the band 
flatness around the (7r,0)-point. Some peculiar features of 
the 2-pole solution (e.g., the inflexion around the (7t,7t)- 
point) can be now clearly interpreted as due to the neces- 
sity of miming the behavior of both Hubbard subbands 
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FIG. 6: Same of Fig. gjbut at U/t = 8. 



and shadow bands by means of only two bands. 

In Fig. [3 we present the doping dependence of the den- 
sity of states and of the corresponding dispersion relation. 
The Fermi level stays pinned to the band flatness located 
at (7r,0)-point within a wide range of hole-doping (0 < 
S < 0.15). For higher doping, it moves towards higher 
energies giving an electron-like Fermi surface centered 
at (0,0)-point. Those characteristics are commonly ob- 
served in many theoretical ana ly S i322i22i2&i2ii2£i2£i which 
take into account nearest-neighbor spin and charge fluc- 
tuations, in agreement with QMC dataiiS^ Therefore, 
we can conclude that in order to reproduce such pecu- 
liar features, it is necessary to take into account high- 
order operators in the basis, as we have done in this 
manuscript. 

In Fig. El we provide a detailed comparison of the band 
structure obtained by the present formulation with QMC 
results. As can be easily seen, we have a good agree- 
ment with QMC data, especially for the low-energy band 
around the Fermi level. We observe shadow bands more 
pronounced than in QMC data. We should recall that the 
present formulation is a pole-approximation and damp- 
ing effects are neglected. Furthermore, three-site terms in 
the equations of motion have been neglected. Therefore, 
there is no diffusion process to weaken two-site correla- 
tion effects. On the other hand, we can expect that with 
hole-doping the antiferromagnetic correlations are weak- 
ened by hole motion and that the shadow bands become 
broader and broader owing to damping effects. Eventu- 
ally, shadow bands are wiped out and we may observe 
traces of them as shoulders of the main Hubbard bands, 
as seen in QMC results. 



1.8 
1.6 
1.4 
1.2 
S 1 



0.2 



1 1 1 1 
~ U/t=8 


i i i i i 








— 5=085" 




- — ~~ ^J5- 


"if 

- / if 




If/ 7 

■l! / 

i/i/ 

i i 
■i i 


- 


■ i 

i , i 






; — i — i — i — | — i — i — i — | — i — i — i — | — i — i — i — | — i — i — r- 



T/t 



FIG. 7: Entropy at U/t — 8 (solid line) (Finite temperature 
Lanczos data (dashed line) from Ref. l4Cfl . Data are shifted 
of 0.2 along the vertical axis for the sake of clarity. 



IV. THERMODYNAMIC QUANTITIES 

Figure 0] presents double occupancy D = (n^ni) = 
(vi 7 h) m comparison with Lanczos^ and QMG^S, results. 
For U/t < 3, it is difficult to get self-consistent solution 
as split bands start merging. Our results show a good 
agreement in a wide range of U jt values. We should point 
out that the 2-pole approximation also provides similar 
agreement^ 

Internal energy E = (H)/N and specific heat C = 
dE/dT per site at U/t = 4 and 8 and n = 0.75 and 0.90 
are reported in Figs. and [S] Data from finite tempera- 
ture Lanczos^ and QMG— are provided for comparison. 
It is worth mentioning that our results for U/t = 12 have 
the same general features that those for U/t = 8, but 
with more pronounced characteristics. As regards the 
internal energy, the agreement with the Lanczos data is 
excellent except for the low temperature region T/t < 0.4 
at U/t > 8. As regards the specific heat, we observe a 
sharp peak around T/t ~ 0.3 and a fairly broad peak in 
the higher temperature region T/t > 1. The two peak 
structure is more pronounced for U/t > 8, but not so 
much for U/t = 4. This tendency is also observed in 
several numerical simulations42^i*^ Usually, the sharp 
peak at lower temperatures and the broad peak at higher 
temperatures are interpreted as consequences of spin and 
charge fluctuations related to the energy scales of J and 
[/, respectively. The main difference between our results 
and numerical ones regards the height of the peak in the 
specific heat around T/t ~ 0.3 that comes from the de- 
crease in the internal energy. This is an indication of 
well established spin ordering which cannot be correctly 
evaluated on a small cluster. Numerical simulation can- 
not describe spin and charge ordering in the case that 
the correlation lengths exceed the cluster size. On the 
other hand, in our formulation, as already discussed in 




FIG. 8: Entropy vs n at U/t = 4, 8 and 12 and T/t = 0.5, 1.0 
and 3.0 (lines) (Finite temperature Lanczos data (markers) 
from Ref. E3). 



Sec. Mil there is no diffusion process to weaken two-site 
correlation effects. Therefore, there is a tendency to have 
too pronounced spin and charge correlations. 

We can discuss this issue in more detail by commenting 
our results on entropy 



S(T,n) = - f 
Jo 



dn 



df 



(57) 



This relation is derived from the thermodynamic 
relations^ S = - (§£) n and \i = (|£) T , which give 

the Maxwell relation (§§) T = - (§f \ ■ Figure re- 
ports the results of entropy in comparison with Lanc- 
zos ones. We present results only for U/t — 8, but 
U/t = 4 and 12 results show the same tendency. Our 
results completely coincide with Lanczos ones in the high- 
temperature region (T/t > 1). In the low-temperature 
region (T/t < 1), our data are lower than Lanczos ones 
indicating a stronger tendency to ordering. 

To better understand the tendency to ordering, we 
have investigated the filling dependence of entropy at 
several temperatures (see Fig. |SJ). For T/t > 1.0, the 
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agreement is extremely good in whole range of filling. 
However, at lower temperatures, our results show a de- 
crease around both quarter and half filling. Usually, a 
decrease of entropy at quarter and half filling is inter- 
preted as an indication of charge and spin ordering, re- 
spectively In a small cluster, it is difficult to investigate 
such ordered states because of the system size. 

V. SUMMARY 

In the present paper, we have carried out an analysis of 
the two-dimensional Hubbard model by means of a 4-pole 
approximation within the Composite Operator Method. 
Density of states and corresponding dispersion relation 
show remarkable characteristics: four-bands structure, 



a quasi-particle peak at the Fermi level, shadow bands 
and band flatness at (71", 0)-point. Quantities such as 
double occupancy, internal energy, specific heat and en- 
tropy have been comprehensively investigated. Our re- 
sults show an excellent agreement with numerical simu- 
lations except for the low energy features related to spin 
and charge ordering. In the present formulation, nearest- 
neighbor-site effects are probably over-estimated whereas 
on a small cluster there are some difficulties to describe 
spin and charge ordering in the case that the correlation 
length exceeds the system size. Probably the right stays 
in the middle: our results provide meaningful informa- 
tion for band structure and thermodynamic quantities, 
but the inclusion of damping effects is necessary for com- 
plete understanding of this system. 



APPENDIX A 



The equations of motion of £ S a(i) and rj SIJ (i) reads as 
d 

' -n- a (i)cf (i)+n_ ff (i)(n_ CT (i)c«(i)) a -n_ CT (i)(c t _ (7 WcaWc^W) Q 
+n^(i)(c a _i(i)c^)c^{i)r cl„(i)c1„(i)tf(i) + c a _l(i)c- a (i)£(i) 
-At I + c -A*M)c- 2 Ai) - d a (i)c„(i)(n„(i)c%(i)r + cl ff (i)<v(f)(4(0^«r(0cS(0) 
] -cl a (i)c„(i)(c?(i)c- a (i)c„(i)r + cl a {i)c«{i)^) - c a J a (i)c a (i)^ a (i) 
-»7-t(iMi)<£<r« - V-Ui)c2(i)c-„(i) + (4 (i) c L <r (i) cS !(i))« Ctr (i) c _ (7 (i) 
{ +(c a J a (i)n a (i)) a c a (i)c^(i) + (c?(i)c a (i)cl a (i)) a c a (i)c_ a (i) 



i-T^Vsvi-i) = (-/" + U)r} srT (i) 

-cl a (i)c%{i)^{i) + c a \(i)c^(i)r,«(i) + cL a (i)c„(i)(n„(i)c%(i)) a 
-c t _ CT (j) Cff (*)( C t(j)c- (T (j) C «(j))«+ C t _ CT («) Cff ( J )( C «t(j) c _ (T ( J ) Ca ( J ))« 

+cl a (i)c%(i)r ] « a (i) - c a J a (i) Cr7 (i) V - a (i) - eUi)c„(i)c a -Ai) - CU^)c-^) 
I -(c«t(;) CCT (;)cl CT (;))« CCT (i) c _ a (i) 



-At I 



where A a (i) and A a (i)B a (i) stand for 



and 



A a (i)B a (i) = ^ayA(i)a ik S(fc) 



respectively. We can isolate the single-site terms as follows 

A a \t) = \A(t) + J2^ ik A(k), 



(A2) 



(A3) 



(A4) 



(A5) 



j#k 
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A a (i)B a (i) = \ (A(i)B(i)) a + J2 ^A(j)a ik B(k). (A6) 

Then, we can isolate the one- and two-site terms in the equations of motion and neglect high-order terms for 
simplicity 



d fl ' 

i— 77^ (*) ~ (-/x + E/)»7 a<r (*) + 4t -^(i) + 

Recalling the relations between Vv(i) and ip a (i) 

U(i) = U(*)-4n(-*V)&(*) 
j? SCT (i) = ?7 SCT («) - A 42 (-iV)»? (T (i), 



we have 



-4i 



- (a(-iV) - A 42 (-zV)) A 3 i(-iV)^(i) - Q - A^-iV^-iV)) »&,(*) 

- (a(-iV) - Ai2(-*V))6 <r (i) + ^ 42 (-«V)% CT (i) 



Then, e BB (k) can be obtained by simple inspection 



where A13 and A24 are defined in the main text. 



(A7) 



(A8) 



d - 

-(A 31 (-iV)a(-iV)+A 31 (-i\7)A 31 (-iV) + a(-iV)A 31 (-iV))Z4i) 
-41 { - Q + A 31 (-iV)a(-iV) + A 31 {-iV)A 42 (-iV) + 2a{-iV)A 42 (-iV^j Vtr (i) 

-(a(-iV)+A 31 (-iV))U(i)-(2a(-iV) + A 31 (-iV))fj s <r(i) ' ) (A9) 



e 00- I M + 4i(a(k) + A 31 (k)) 4t(2a(k) + A 3 i(k)) \ 

SBl ; I 4i(a(k)- A 42 (k)) -,j + £/-4L4 42 (k) j' 1 j 



APPENDIX B 



7 3 3(i,j), ^34(1, j), and 1*4 (i,j) reads as 



= 61 



WU) = ({Ui). 

/ (d CT (i)e CT (o«- CT (i)> - (^t(0»?-«r(i)n-«r0')> + (Ct(j)e-.o>-,») \ 

( (n- a {i)n- a (j)) - (n_ (7 (i)n_ (7 (fc)n_ <r (j)> - (cL CT (i)cv(i)4 (fc)c- CT (fc)n_ CT (i)} \ 
+<c IT «c_ CT (i)c t _ CT (fc)4 (fc)n_ ff (j)) - (n_ a {i)cl a {k)c a {k)cl (j)c_ CT (j)} 
+{J_ IJ (i)c a {i)c] a {j)c- a {j)) - (cL„(i)c a (i)n„(k)cUj)c-v(j)) 



(Bl) 
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Wi,j) = <{£«r(M), 4Ai,t)}) 



= 6; 



(B2) 



^44a(i,j) = ({r/ sff (M), »?L(j,t)}> 

/ (£l(i)^(i)) (£\{i)Z a -Ai)n a {i)) 

V -(V-l(i)v2(i)cl (i)c-„(i)) (v?(i)£Ui)cAi)c-Ai)) + {clMWtfWP,®) 
( -(Ci(i)Z-a{i)n- a (j)) + ( V l a {i) V %{i)n. a (j)) + (eUj)V-*(j)n-a(i)) \ 

( (n_ CT (i)n_ CT (fc)n_ CT (i)) + (cL CT (i) ClT (i)4 (k)c_ a (k)n_„(j)) \ 
+(cl a (i)c a (i)n a (k)cUj)c-a(j)) + (n- a (i) c -A k ) c A k ) c i (j)c- a {j)) 
-(c a ( J )c_ ff (i)cL (T (fc)ct(fc)n_ (T (j)) + ( C t _ (T (j)4(j>- ff (fc)ca(fc)n- ff ( i )) 
V +(cL CT (j)4 (j)c CT (i)c_ CT (i)> - {ciMH(j)nAkM)c-«W) ) 



+a i} 



(B3) 



aikCKkj 

k 



They contain three-site correlation functions which cannot be directly evaluated in terms of the propagators under 
analysis. Therefore, we have decided to decouple them in terms of two-site correlation functions. For example, 



< C U*)4(*)C«^(*)> 

^hcl a (i)ct(i) (e ff Wr ? _ CT (z)) a )+^a ij « ik ( C t_ (T ( i )^(fc)}(4(zK CT (. 7 )). 



(B4) 



Within this procedure, those terms reducible to one- and two-site correlation functions have been exactly evaluated. 
Because of the translational symmetry, we have 



f 33CT (k) ~ \ «n_„) -p a ) + fa-U-«,> - 



+a(k) 



-2A ff |(n_ CT ) + 1 (1 - («_„))} - 2<£ (7 <£ <r X££0 + ^Otfl^-*) 
+2(^ ffC « J^c-J - \<rf_ a <?_ a )(r{L a n-*) 2(4^)^4) + \(cl a c%)(tU a ) 

+2( v t^)( v a _U a _ a ) l(vt^)(vUv-a) 2( v tO(CU-,) + livteKZUt-*) 



(B5) 



r/3 2 (k) 



(n^)(n^) (1 - <n_ ff )) + 2(n_ (T )(4_ CT c« <T ) ((cl^J + 2(44) 
-(cL.^Xct.^Xl - 2 (n_ CT )) - (c t _ ffC ^)(44 1 )(l - 2<n CT » 
-2((44) + (c t _ ffC ^)) 2 (c t _ (TC ^) 

(n_ CT )(n_ CT ) (1 - (n_ ff » + 2(n_ <7 )(cL CT c« <7 ) ((d^J + 2(44)) 
-(cL^Xcl^Xl - <«_„» - {cl a cll){cU^)(l 2(n a )) 



((44) + (c^j) (cl.c* > 
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Wk) ~ ( V a _U« a ) (1 - 2<n_ CT » + (l - J) (2tf-M(<L*g.*) + {4<%)(ctM) 
+a(k) f . v . ' 

mo bus)) + 4 - (^^-.)) - bus)) 



(B6) 



/44.(k)^rU- CT > + ^^-.) 



+ 1 



+a(k) 



+ i /3 2 (k) 



-<cl ff c* X«t AX"-*) + (c t - ff c^)( c t c ft)(l - 2<n„» ' 
+2((ct c «) + ( c t_ CTC « ff )) 2 ( c t_ c ft) 

<n_„Xn-*X»-<r> - 2(n_ (T )(c t _ (T ^ CT ) ((cl^) + 2( C t c «)' 
-{cUc* Xct^c* Xn.,) + (ct.c* Xctc^)(l - 2(n CT )) ' 



(B7) 



where 



with 



a(k) = — {cos k x a x + cos k y a y } 

/?i(k) = cos k x a x cos k y a y 

/?2(k) = - {cos 2k x a x + cos2k y a y } 



(B8) 



a 2 (k) = i + i/3 1 (k) + i/? 2 (k). 



(B9) 



Parameters p<j and A CT arc defined in the text. 

{A^B^ 1 (i)) and (A(i)B^ 2 {i)) contain next-nearest-neighbor-site correlation functions along the diagonal and the 
main directions, respectively. If we assume that those correlation functions have same value, that is, we use the 
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spherical approximation, we can simplify the momentum dependence of I 



-2A CT |(n_ CT ) + i (1 - (n_ ff ))| - 2<£l £r c« ff XetO + ^O^-*) 



+a(k) 



+/3(k) 



+2{r ? t c «)(^t CTC a (r) _ i^)^^) - 2(r ? t c «)(^ t CT tf a ) + ~<^c?)(d CT £- CT ) 
(n_ CT )(n_ CT ) (1 - + 2(^)^0 f (cl CT c« CT ) + 2< C t c «0 



-2(( C t c «) + ( c t_ CTC « CT )) 2 ( c t_ CTC ^) 



(BIO) 



with 



+a(k) 



+/3(k) 



2 1 



1 



1 



(n- ff )(n_ CT )(n_ ff ) - 2{n^)(cl a c a _ a ) {{cL a c« a ) + 2( C T C «) 
-(cU^XcU^Xn.,) + (cU^XctcgXl - 2(n CT )) 

+2((4 c «) + { c L (tC ^)) 2 ( c L ctC ^} 



/3(k) = i/3 1 (k) + i/3 2 (k) = a 2 (k)-i. 



(Bll) 



(B12) 



We have checked that this assumption doesn't produce any difference as regards the results reported in the present 
paper. We have used this approximation in order to simplif y the momentum integration. 
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